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Mott Transitions and d-Wave Superconductivity in Half-Filled-Band Hubbard Model 
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Mechanisms of Mott transitions and d^2_y2-wsLve superconductivity (SC) are studied in the 
half-filled-band Hubbard model on square lattices with a diagonal hopping term {t'), using an 
optimization (or correlated) variational Monte Carlo method. In the trial wave functions, a 
doublon-holon binding effect is introduced in addition to the onsite Gutzwiller projection. We 
mainly treat a d-wave singlet state and a projected Fermi sea. In both wave functions, first- 
order Mott transitions without direct relevance to magnetic orders occur at U ~ Uc, which is 
approximately the bandwidth, for arbitrary t' /t. These transitions originate in the binding or 
unbinding of a doublon to a holon. d-wave SC appears in a narrow range immediately below 
Uc- The robust d-wave superconducting correlation is necessarily accompanied by enhanced 
antiferromagnetic correlation; the strength of SC decreases, as t' /t increases. 

KEYWORDS: Mott transition, superconductivity, condensation energy, antiferromagnetic correlation, spin 
gap, high-Tc cuprate, Hubbard model, variational Monte Carlo method, frustration 



1. Introduction 

In connection with the superconductor-insulator tran- 
sitions in organic compounds k-(BEDT-TTF)X,^ the 
half-filled-band Hubbard model on anisotropic triangular 
lattices [Fig. l(b)]^ has been intensively studied. In ad- 
dition, a recent experimental study'^ has reported that a 
series of film samples of nondoped high-Tc cuprates (par- 
ent materials of electron-doped systems) do not become 
antiferromagnetic (AF) insulators but exhibit metallic 
properties including superconductivity (SC) at 21 K. 
Thus, it is important to grasp the mechanisms of Mott 
transitions and SC, if any, in half-filled-band Hubbard 
models with frustration. Because such phenomena arise 
at intermediate correlation strength {U ~ W; W be- 
ing the bandwidth), one has to use a method that can 
reliably treat both strongly correlated and weakly corre- 
lated regimes. As one of such methods, the optimization 
(or correlated) variational Monte Carlo (VMC) method* 
has rapidly progressed in recent years for the study of 
ground-state properties.^ 

Using this method, the present authors have recently 
studied the Hubbard model on the anisotropic triangu- 
lar lattice [Fig. l(b)].^ With a projected d-wave singlet 
state, it is found that (1) a conductive-to-nonmagnetic- 
insulator (Mott) transition occurs at U — Uc, which is 
somewhat smaller than W, and is caused by the binding 
(and unbinding) of a doublon to a holon; (2) c?2,2_y2-wave 
SC appears in the vicinity of both the Mott transition 
and the AF phase, and develops together with the short- 
range AF correlation (or fluctuation). 

In weakly frustrated cases, the two lattices in Figs. 1(a) 
and 1(b) have a common characteristic wave number 
G — (tt, tt) in spin correlation. In strongly frustrated 
cases, however, the characteristic wave numbers become 
different; for instance, 120-degree-structure spin correla- 
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tion is probably dominant in (b) for t' ~ t,^ whereas 
the collinear structure is favored in (a). In this paper, 
we carry out similar detailed calculations for the lattice, 
often treated in the context of high-Tc cuprates [t-t'-U 
model) [Fig. 1 (a)] , and reveal whether or not the above 
mechanisms of the Mott transition and of SC also work 
here. In addition to the d-wave singlet state, we study 
Mott transitions in the projected Fermi sea. 
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Fig. 1. (Color online) Lattice structure and hopping integrals t 
and t' , (a) studied in this work, and (b) often used for k-BEDT- 
TTF salts. Lattice sites are denoted by dots. 



The organization of this paper is as follows: In §2, we 
introduce the model and method we use. In §3 and §4, 
we discuss Mott transitions in the c?-wave singlet state 
and in the projected Fermi sea, respectively. Section 5 is 
assigned to the stability and properties of the d-wave su- 
perconducting (SC) state. In §6, we construct a ground- 
state phase diagram based on the present VMC calcu- 
lations, and address important problems with respect to 
antiferromagnetism (AF). In §7, we briefly summarize 
the main results and discuss the subject further. 

Part of the present results have been reported in a 
previous letter.^ 

2. Model and Method 

In §2.1, the model we study is introduced, and related 
studies are summarized. In §2.2, we briefly review the 
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background of variational wave functions for the Hub- 
bard model in the research on the Mott transition. In 
§2.3, we give an account of the wave functions used in 
this paper and the conditions of VMC calculations. 

2. 1 Hubbard model on frustrated square lattices 

In this paper, we study the Hubbard modeP^^^ on a 
square lattice with diagonal transfer t' [Fig. 1(a)], 

U = Hicin+Hu = ^£(k)c[^CkCT + U ^fij^riji, (1) 



kCT 



e(k) — — 2i(cos + cos ky) — 4t' cos k^ cos ky 



(2) 



with ?7, t > 0. Equation (1) has been often used as a sim- 
ple model that captures the essence of cuprates.^^ Here, 
we concentrate on the half- filled band {n — N^/Ns = 1; 
Nc- electron number and Ng'. site number) and consider 
doped cases in a forthcoming publication. We exclusively 
treat the case of t' < 0, because the behavior for t'{> 0) 
is identical to that for —t', owing to the particle-hole 
symmetry at n — 1. Note that the negative sign of 
t' /t is in agreement with the hole-doped case of high- 
Tc cuprates. When \t' /t\ is varied, the features of the 
bare band, eq. (2), abruptly change at \t'/t\ = 0.5, at 
which low-lying energy levels become strongly degener- 
ate; the van Hove singularity points depart from (tt, 0) 
and (0, tt) for \t' /t\ > 0.5. Furthermore, the plausible val- 
ues of high-Tc cuprates are considered to be \t' /t\ = 0.1- 
O.S.^'^ Thus, we restrict the range of frustration strength 
to < 1^7*1 < 0.5 in this paper. 

Despite the importance of the model, reliable knowl- 
edge is limited, particularly in the intermediate and 
strong coupling regimes. For the pure square lattice 
{t' — 0), it is believed that the ground state is insulating 
with a long-range AF order for U > 0, owing to the com- 
plete nesting condition.^"* For frustrated cases {t' ^ 0), 
however, it is urgently required to clarify the properties 
of the conductor-insulator transition. For SC at half 
filling, although the anisotropic triangular lattices have 
often been studied, studies on the present lattice are 
rare to our knowledge. 

2. 2 Historical background of wave functions 

As a many-body trial wave function, the Jastrow 
type^" is useful and has been often applied: ^' = 7^$, 
where V denotes a many-body correlation (Jastrow) fac- 
tor composed of projection operators, and $ is a one- 
body wave function usually given by a Slater determi- 
nant. For the Hubbard model, more than four decades 
ago, Gutzwiller introduced the celebrated onsite projec- 
tion, 
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- n [1 - (1 - 9)nnnji] 



(3) 



which has primary importance for arbitrary parameters 
in the Hubbard model. Although the Gutzwiller wave 
function (GWF), ^-^^ = Vg^f {"^f- Fermi sea), looks 
simple, it is generally difficult to accurately calculate 
expectation values using it. Hence, a mean-field-type 
approximation [now called a Gutzwiller approximation 
(GA)] was introduced by Gutzwiller himself,^^ and was 



used and extended by many researchers for the following 
two decades. However, variation theory loses its vari- 
ous advantages when additional approximations such as 
GA are applied, and consequently it becomes difficult 
to improve the wave function. To break this deadlock, 
VMC methods^'^ have been applied to this problem;^^' 
thereby and by subsequent exact analytic treatment in 
one dimension, the precise behavior of ^fg^ was clari- 
fied for the first time. Although Vq is indispensable for 
treating the Hubbard model, its independent use leads 
to the following physically unsatisfactory results: (1) The 
momentum distribution function n(k) tends to be an in- 
creasing function of |k|, (2) 2kp anomalies in the spin 
[charge density] structure factor 5(q) [A^(q)] cannot be 
properly represented, and (3) a Mott transition cannot 
be described, in addition to the problem of a considerably 
high variational energy. 
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Fig. 2. (Color online) Weight assignment of Jastrow factor V = 
VqT'g depending on local electron configuration. Each circle in- 
dicates a site. A solid (open) circle indicates a doublon (holon). 
Thin arrows denote virtual hopping processes in the strong- 
coupling expansion, (a) Configuration with no doublon; a basis 
for U/t — > oo. (b) A doublon sits at a nearest neighbor of a holon; 
a virtual state in the second order of strong-coupling expansion 
in t/U. (c) A doublon sits at a diagonal neighbor of a holon; a 
virtual state in the second (fourth) order in t' /U (t/U). (d) A 
doublon sits at a farther site from a holon; a higher-order virtual 
state. For the case of eq. (4), we set /i' = 0. 



Although the electron-electron interaction in the Hub- 
bard model is limited to within a single site, its effect 
reaches distant sites. Therefore, to overcome the above 
shortcomings of Vq, one needs to add intersite correla- 
tion (long-range Jastrow) factors. In cases of low electron 
density, distance-dependent long-range Jastrow factors 
are useful. On the other hand, at half filling, the short- 
range part of the Jastrow factor is predominant owing to 
the screening effect. Castellani et al."^^ derived an effec- 
tive Hamiltonian of the Hubbard model, taking account 
of both spin and charge degrees of freedom. In their ef- 
fective Hamiltonian, an exchange term between a dou- 
bly occupied site (doublon) and an empty site (holon) 
appears, indicating that a doublon- holon correlation is 
inherent in the Hubbard model. Kaplan et al?'^ actually 
showed by studying one-dimensional (ID) small clusters 
that the binding of a doublon to a holon is important for 
large U/t to reduce the energy at half filling. Using ex- 
act diagonalization, Yokoyama and Shiba^^ studied the 
ground-state wave function of Hubbard rings at half fill- 
ing, at which the ground state is known to be insulating 
for U/t > 0.^° They found that, for large U/t, the magni- 
tudes of the coefficients of bases with one doublon (and 
one holon) decrease exponentially as a function of the 
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distance between doublon and holon. This means that a 
doublon is bound to a holon within the decay distance 
in an insulating state. 

A simplified wave function that reflects the above 
arguments is written as = VqVq^-^'^''^^ Here, 
the doublon-holon binding factor Vq is limited to the 
nearest-neighbor part: 



^Q=n(i-'"^n' 



(4) 



2.3 Wave functions and VMC conditions 

In this paper, we continue to study the Mott transition 
induced by Vq. Considering the lattice structure shown 
in Fig. 1(a), we introduce into Vq the effect of doublon- 
holon binding between diagonal-neighbor sites /i' [r = 
(±a;, ±2/)], in addition to that of the nearest neighbors, 
fj,[r= {±x,0) and (0,±y)]: 



rQ = l[{l-f^Qj){l-f^'Qj'), 



(7) 



(5) Q 



r(r') 



where di = Ui-^Uii, Cj = (1 — nii){l — Uii), and t varies 
over all the nearest neighbors. A variational parameter 
M (0 ^ M ^ 1) controls the strength of doublon-holon 
binding in the nearest-neighbor sites; as /x increases, a 
doublon tends to adhere to a holon, and for — > 1, a 
doublon cannot leave a holon, as shown in Fig. 2. 

This Jastrow factor, VqVq, can also be derived nat- 
urally from the strong-coupling expansion. It is known 
that the GWF with 9 = is an extremely good trial 
state for the ID Heisenberg model, ^^'^^ and that the 
Gutzwiller-type functions, 'Pg^af and 'Pg^bcs (^af: 
Hartrcc-Fock-typc AF state: "I^bcs: (i,,2_j,2-wave BCS 
state) , yield quantitatively reasonable results for the 2D 
t-J model.^^"^^ These favorable properties of the GWF 
for strong-coupling models can be applied to the Hub- 
bard model by considering a canonical transformation, 
Ut-j ~ e'-^Wnube"'^ {t/U 0),^^ and similarly. 



(^G|H,-,yhI^G) _ ('^G':'lHuubk""^I^G) 

^ (*Ge*'^|e~*'^*G) 



Thus, an improved wave function for the Hubbard model 
is given by applying the strong-coupling expansion, e~^^ , 
to a Gutzwiller-type function (= Vq^). Considering 
virtual hopping processes in this expansion, as shown 
in Fig. 2, one easily notices that the first-order terms of 
this expansion roughly correspond to Vq^q. In addition, 
it was found that a more direct form of e~*'^4'G yields 
improved results"^^ similar to those oiVq^G, mentioned 
in the following sections. 

In the early VMC study of the projected Fermi sea, 
*FS ^ PqPg$f, for the ID and 2D square lattices, 
Yokoyama and Shiba^'^ concluded that corrects the 
shortcomings (1) and (2), mentioned above, of the GWF, 
but a Mott transition does not arise, even in Sub- 
sequently, Millis and Coppersmith^^ also concluded, by 
calculating the zero-frequency part of the optical con- 
ductivity using a VMC technique, that a Mott transi- 
tion cannot be described in terms of this type of wave 
function. However, as we have repeatedly explained in 
previous papers,^' these early studies were not suffi- 
ciently thorough or careful to arrive at the correct con- 
clusion that the doublon-holon binding factor Vq is es- 
sential for describing a Mott transition. Actually, the ex- 
istence of Mott transitions has been confirmed using Vq 
for various systems, "^^ such as the square lattice,^' the 
anisotropic triangular lattice,^ the kagome lattice, **° the 
checker-board lattice,^^ and a degenerate Hubbard model 
on the square lattice. 



n 

r{r') 



ei+T{T')) + ej(l - di+r{r'))\ , (8) 



in which t (r') varies over all the adjacent sites in the 
bond directions of t {t'). The weight assignment of the 
correlation factors V = VqVg is explained in Fig. 2. For 
the pure square lattice {t' = 0), wc use cq. (4) instead 
of eq. (7) for simplicity, because we can confirm that the 
effect of jj.' is negligible, even quantitatively, for t' = 0. 
The wave function we deal with in this paper is 



(9) 



As a one-body part we primarily study a fixed- 
density BCS state:^3 



|0), 



Ok 



£k-C + V(£k-C)'+A^' 



(10) 



(11) 



where C is a variational parameter that is reduced to 
the chemical potential for U/t ^ 0. Since we know that, 
at half filling, the simple dx2_y2 wave is the most stable 
among the various gap shapes,^^'^^'^^ here we exclusively 
consider the d-wave gap: 



Ad (cos kx - cos ky ). 



(12) 



Note that although the variational parameter indi- 
cates the magnitude of the d-wave gap, a state {^q = 
V^d) with finite A^ does not necessarily mean a SC 
state. ''^ We fix the value of t' in ek in the wave func- 
tion [eq. (11)] at the same value as that in the Hamilto- 
nian, because the renormalization of £k^^ is not signifi- 
cant when the system is conductive,^ and ( compensates 
this effect to some extent in the insulating regime, as we 
will see later. For A^ = 0, $d is reduced to $f, which is 
explained next. 

As a reference state, we also study the Fermi sea, 

n 4jo). (13) 

k<kF,a 

A complication is that ^fg^ (= VqVg^fs) does not 
merely represent a normal state; also undergoes a 
Mott transition, as we found for the attractive Hubbard 
model. The results of the attractive model for the sym- 
metric case {t' = and n = 1) can be exactly mapped to 
those of the repulsive model through a canonical trans- 
formation.'*^ As an extension in this paper, we study the 
properties of \E'q^ for asymmetric cases {t' ^ 0). 
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For the comparison made in §5 and §6, we consider an 
ordinary mean-field solution $af^^ for a one-body AF 
state with a long-range order. In the trial AF state, 



G'i'AF, 



(14) 

the AF gap Aaf is optimized as a parameter, but t' in 
£k is fixed at the model value, similarly to that in ^'q. 

Because our trial functions have at most five parame- 
ters to be optimized [g, /i, /i', A^(af), Q], we have used a 
simple version of optimization VMC methods,"*^ namely 
a line minimization of one parameter with the others 
fixed. In one round of iteration, every parameter is op- 
timized once. In most cases of optimization in this study, 
the parameters converge within 2 or 3 rounds, after which 
we continue the optimization process for another 15-20 
rounds. The optimized values of the parameters and en- 
ergy are determined by averaging the results of these 
rounds after convergence. Because each optimization pro- 
cedure is carried out with typically 2.5 x 10^ {L = 10- 
14) samples, preserving the acceptance ratio of 0.5, our 
data are practically the averages of 3-5 million samples. 
Thereby, the accuracy in the total energy is markedly in- 
creased, typically to the order of 10~'^t. Because the con- 
vergence of optimization becomes very slow near phase 
transitions, particularly continuous transitions, we car- 
ried out longer iterations (~ 50 rounds) in such cases. 
With the optimized parameters thus determined, phys- 
ical quantities are calculated in different VMC simula- 
tions with 2.5 X 10^ samples. We used lattices oi L x L 
sites (L — 6-18, mainly 10-14) with periodic-antiperiodic 
boundary conditions. Because the closed-shell condition 
cannot be satisfied in a wide range of asymmetric model 
parameters {t' ^ 0) at half filling, we are often obliged 
to calculate with open shells. 

Finally, we mention finite-size analysis in this study. In 
the symmetric case (t' = 0), the system-size dependence 
of various quantities is often monotonic, because the k- 
point structure included in the Fermi surface becomes 
systematic with increasing L. In contrast, in asymmetric 
cases, the k-point structure is unique to each system size 
L and frustration strength t'\ therefore, the system-size 
dependence becomes irregular, and the t' jt dependence 
is not smooth. 

3. Mott Transitions in d-Wave State 

In this section, we study Mott transitions arising in 
the d-wave singlet state VPg. In §3.1, we show that a 
first-order transition occurs by observing hysteresis in 
the U jt dependence of total energy. In §3.2, we identify 
this critical behavior as a Mott transition without rele- 
vance to magnetic orders by studying various quantities. 
In §3.3, we consider the properties of this transition with 
reference to other studies. 

3.1 Total energy and hysteresis behavior 

First, let us consider the behavior of the optimized 
total energy per site, E. In the inset of Fig. 3, E/t for 
t' /t = is plotted for a wide range oiU/t and four sys- 
tem sizes [L). On this scale, L dependence is impercep- 
tible. One readily notices a cusp at U /t ~ 6.5, where the 
data points of various L are concentrated. In the main 



panel of Fig. 3, this cusp is magnified. We have opti- 
mized the wave function successively, from both weak- 
and strong-correlation sides toward this cusp. For the 
system with L — 10, the optimized energies from the 
two sides coincide, and E/t becomes a smooth function 
of U /t. On the other hand, for larger systems {L > 12 
in Fig. 3), the optimized values in the weak-correlation 
side are not smoothly connected to those in the strong- 
correlation side, and vice versa. At the cusp, all op- 
timized variational parameters have discontinuities, as 
shown in Figs. 4(a)-4(e). Such hysteresis and discontinu- 
ities indicate that a first-order phase transition occurs 
at the cusp point ([/ = Uc). In the following subsec- 
tion, we identify this transition as a Mott (conductor-to- 
nonmagnetic-insulator) transition. 
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Fig. 3. (Color online) Total energies of d-wave state for pure 
square lattice (t' = 0) for four system sizes near critical points 
(Uc), indicated by arrows. Hysteresis is observed for L = 12-16; 
the local minima for the conductive (insulating) sides are denoted 
by open (solid) symbols. For L = 10, hysteresis is not observed 
using the present VMC calculations, and the parameters vary 
continuously. The critical values are Uc/t = 6.54, 6.61 and 6.69 
for L = 12, 14 and 16, respectively. Inset: The behavior of E/t 
over a wider range of U/t. 




Fig. 5. (Color online) Total energies for four values of t'/t near 
critical points indicated by arrows. The system size is fixed at 
L = 12. Hysteresis is observed for each case. The local minima 
for the conductive (insulating) sides are shown by open (solid) 
symbols. The critical values are Uc/t = 6.54, 6.62, 6.70 and 7.00 
for t'/t = 0, -0.3, -0.4 and -0.5, respectively. 





Fig. 4. (Color online) (a)-(c) Optimized variational parameters (at the global minima of E/t) for d-wave state near critical values Uc 
as function of U/t: (a) Onsite correlation (Gutzwiller) factor, (b) d-wave gap parameter, (c) Doublon-holon binding factor between 
nearest-neighbor sites, (d) The same factor between diagonal-neighbor sites, (e) Chemical potential. For t' = 0, C is always zero owing 
to the band symmetry, (f) Total energies for four values of t' /t. The symbols and abscissa scales are common to all the panels. Data 
for four values of t' /t, and for three or four system sizes (L X L) are compared. 



Next, we study the L dependence of the critical point. 
We first note that it is essential to check the system- 
size dependence when we consider critical phenomena, 
as we have learned from the distinct behavior between 
i = 10 and L > 12 in Fig. 3. As the system size increases 
[L > 12), the critical point shifts to a larger value of U/t 
(Fig. 3). This is partly because a large system size is more 
advantageous to conductive states, which have longer 
correlation lengths. For t' /t = 0, because the system-size 
dependence of E/t is fitted well with quadratic curves of 
1/Ns, the critical value for L = oo can be estimated as 
Uc/t = 7.0 ±0.1 using the method of least squares. 



Here, we consider the t' /t dependence. As shown in 
Fig. 4(f), the total energy E/t also has a cusp for finite 
t'/t. The behavior of E/t near the cusps is magnified in 
Fig. 5 for four values oit'/t and for L = 12. For U < Uc, 
E/t is significantly reduced as \t'/t\ increases; this be- 
havior is common to the noninteracting case, in which 
the gain in Et' exceeds the loss in Et [Et {Et> ) : the con- 
tribution from the t [t') term to (Tikin)]- In contrast, for 
U > Uc, E/t changes only very slightly; we will return 
to this point in §3.3. Consequently, as \t' /t\ increases, 
the transition point shifts to a larger value of U/t, par- 
ticularly rapidly for large t'/t, although the bandwidth 
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remains constant, 8t, for < t' /t < 0.5. This t' /t depen- 
dence of Uc/t can be verified from the behavior of the 
parameters [Figs. 4(a)-4(e)]. 

In Fig. 6, we show the system-size dependence oi E/t 
near Uc/t for strongly frustrated cases {t' /t = —0.4 and 
-0.5). Here, even the system as small as L = 10 exhibits 
clear hysteresis. The critical value Uc/t tends to increases 
monotonically as L increases, although the extrapolation 
of Uc/t to L = cx) is difficult using the present data, 
because of the nonmonotonic system-size dependence, as 
mentioned. However, we predict that Uc/t for L = oo is 
only slightly larger than those for finite L, because the 
rate of increase of Uc/t with respect to L is similar to 
that for t' /t = 0. 

3. 2 Confirmation of Mott transition 

In this subsection, we confirm that the above transition 
is a Mott transition by studying various quantities. 

First, we consider the doublon-holon binding factor 
/X, which is a good indicator of the Mott transition. In 



Fig. 4(c), the optimized value of /i is plotted as a func- 
tion of U/t. For each t' /t and L, except for t' — and 
L = 10, has a discontinuity and suddenly increases 
at J7 = C/c- Note that for [/ < Uc, primarily de- 
pends on t'/t and slightly on L in relation to SC, as 
mentioned later, whereas for U > Uc, fJ- significantly in- 
creases as L increases, but is almost independent oit' /t. 
Thus, in the strong-correlation side U > Uc, fJ- has a value 
close to 1,^^ which means that almost all doublons and 
holons are bound within nearest-neighbor sites. This sit- 
uation is shown in the snapshots taken in the VMC sam- 
pling process (Fig. 7). On the weak-correlation side of Uc 
[Fig. 7(a)], where has a relatively small value, doublons 
(negative charge carriers) are often isolated from holons 
(positive charge carriers). Thus, charge can move sub- 
stantially, and the system is considered conductive. On 
the other hand, for U > Uc [Fig. 7(b)], each doublon is, 
in most cases, paired with a holon, in addition to there 
being a decrease in the carrier number. It follows that 
free charged particles rarely exist. 




for three system sizes. Hysteresis is observed for every case; the pjg g (^ofor onhne) Density of doublons (doubly occupied sites) 

conductive and insulating cases are denoted by the same sym- obtained with as a function of U/t for four values of t'/t. For 

bols. The critical values are Uc/t = 6.64, 6.70 and 6.78 for ^, ^^^-^^^^ ^^^^^^ ^i^g^ ^ ^Q_^g) simultaneously 

L = 10, 12 and 14, respectively, for t'/t = -0.4, and similarly, plotted. The symbols are common to Fig. 4. For t'/t = 0, the 

Uc/t = 6.93, 7.00 and 7.10 for t'/t = -0.5. horizontal axis is shifted by 0.5 (see upper axis) for clarity. In 



the inset, the same curves are shown over a wider range. 
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(a) U<Uc 








(b) U>Uc 









Fig. 7. (Color online) Snapshots (typical electron configurations) 
taken in VMC sampling process of d- wave state for t'/t = and 
L = 16 (Uc/t = 6.69). (a) U/t = 6.25 (p. = 0.30), slightly less 
than Uc, and (b) U/t = 7.50 (^t = 0.85), slightly greater than 
Uc- Closed and open squares, upward and downward triangles 
denote doublons, holons, f- and J,-spins, respectively. 



Second, we consider the doublon density, 

]^l]'''T"a = ^> (15) 

^ i 

which is regarded as the order parameter of Mott tran- 
sitions,^^' by analogy with the particle density in 
gas-liquid transitions. In eq. (15), Eu = {Hint)/Ns. In 
the inset of Fig. 8, d is plotted for four values of t'/t and 
for various system sizes. As U/t increases, d decreases lin- 
early from the noninteracting value 0.25, but a,t U — U^ 
it suddenly drops to a considerably smaller value, and 
then decreases slowly for U > C/c- In the main panel of 
Fig. 8, the vicinity of the critical point is magnified; the 
discontinuity of d at Uc is clear for each case. This abrupt 
decrease in doublon density can actually be verified in 
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Fig. 7, where the number of doublons is 24 {d — 0.094) 
for U/t = 6.25, and 9 {d = 0.035) for U/t = 7.5, which 
is consistent with the values in Fig. 8. Thus, the discon- 
tinuity of the doublon density at Uc is similar to that of 
the change in particle density in gas- liquid transitions. 




Fig. 9. (Color online) Momentum distribution function of the d- 
wave state for various values of U /t for t' /t = 0. The half-closed 
symbols and star denote the data for U > Uc- The arrow on 
the upper axis indicates the position of the Fermi surface in the 
node-of-gap (F-M) direction. 



Third, the behavior of the momentum distribution 
function. 



Cko 



(16) 



at the Fermi surface is another good indicator of a Mott 
transition. In Fig. 9, we show the U/t dependence of 
n(k) measured with the optimized parameters along the 
path r(0,0)-X(7r, 0)-M(7r, 7r)-r in the Brillouin zone for 
t' /t = 0. Because the present trial state is a projected 
d-wave, there is a node in the gap function in the F- 
M direction, and n(k) has a discontinuity at the Fermi 
surface, kp, in this direction if the system is metallic or 
SC. As shown in Fig. 9, the discontinuity at kp is clear 
for U < Uc, whereas, aX U — Uc, the behavior of n(k) 
abruptly changes, and becomes a smooth function for 
U > Uc also in the F-M direction; that is, the Fermi sur- 
face disappears. Thus, metallic properties are abruptly 
lost &t U = Uc, even in the nodal direction of the dx^^y2 
wave. 

To consider this behavior quantitatively, we em- 
ploy the quasi-particle renormalization factor Z , which 
roughly corresponds to the inverse of effective mass, un- 
less the k-dependent renormalization of self energy is se- 
vere. We estimated Z from the magnitude of the jump in 
?T,(k) at k = kp in the nodal direction, and plotted it in 
Fig. 10 for four values of t' /t.^^ For all the values of t' /t. 



Z decreases slowly for U < Uc, whereas a.t U — Uc, Z 
suddenly vanishes with a sizable discontinuity, reflecting 
the first-order character of the transition. The system- 
size dependence of Z is very small, except for minor dif- 
ferences near Uc- The behavior of Z strongly suggests 
that the effective electron mass diverges for U > Uc- 




U/t 

Fig. 10. (Color online) Quasi-particle renormalization factor Z of 
d-wave singlet states, estimated from discontinuities of n(k) in 
node-of-gap direction. Data for four values of t'/t arc plotted as 
a function of U /t- 




Fig. 11. (Color online) U/t dependence of charge structure factor 
for d-wave state (L = 14) for (a) t'/t = -0.25 and (b) t'/t = 
—0.4. The two panels are in the same scale. The symbols in both 
panels are common to Fig. 9 except for the values oiU/t specified 
here. 
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Finally, let us consider the charge structure factor, 

^(q) = ^5]e^''-(^-^^) {n,n,)-n\ (17) 

with rii = + n^i . Within variation theory, it is known 
that -/V(q) cx |q| for |q| ^ if the state does not have a 
gap in the charge degree of freedom, whereas -/V(q) oc q^ 
if a charge gap opens. In Fig. 11, we show the U/t depen- 
dence of A^(q) for t' /t = -0.25 and -0.4. For U < Uc, 
the behavior of N{q) near the F point seems linear in |q| 
for both values of t' /t. The behavior of iV(q) abruptly 
changes at the critical point, and seemingly becomes 
quadratic in |q| for U > Uc- Although it is not easy 
to definitely determine the power of 7V(q) for the small 
systems used here, we find that the size dependence is dif- 
ferent for U < Uc and for U > Uc- As shown in Fig. 12, 
the behavior of N{q) near the F point for U < Uc ap- 
proaches the analytic curve of U/t = (i = oo) as L 
increases. Conversely, for U > Uc, the slope of A^(q) for 
|q| ^ becomes smaller as L increases, suggesting the 
quadratic behavior of N{q). It follows that ^'g is gapless 
in the charge sector and conductive for U < Uc, but a 
charge gap probably opens for U > Uc and becomes 
insulating. 

Because of the behavior of all the quantities discussed 
above, it is appropriate to judge that a first-order Mott 
transition occurs in the d-wave singlet state at U — Uc 
for the arbitrary values of t' /t considered. 

3.3 Properties of Mott transitions 

Most of the properties of the Mott transition in in 
the present model, eq. (1), are shared with those stud- 
ied in the preceding report for the anisotropic triangular 
lattice.^ In such cases, we avoid repeating detailed expla- 
nations, and only give a brief summary. 

(1) In contrast to the behavior in Brinkman-Rice the- 
ory,^^ in which electrons cease moving and doublons com- 
pletely vanish in the insulating regime, the present VMC 
result exhibits energy reduction broadly proportional to 
-t^/U - J/4) for U > Uc, as shown in Fig. 13. Thus, 
as we argued in a previous letter,^ the results of strong- 
coupling theories (i- J-type models) are qualitatively use- 
ful for U > Uc, the value of Uc is roughly equal to the 




F X q M 



Fig. 12. (Color online) System-size dependence of charge struc- 
ture factor in d-wave state (f /t = 0) near P point for U = 6.25t 
(< C/c) and U = 7.0t (> Uc). 



bandwidth. 

(2) As the frustration t' /t increases, the character of 
the first-order phase transition becomes notable. For ex- 
ample, as t' /t increases, the hysteresis in E/t is observed 
in smaller systems, and the magnitude of discontinuity 
at U — Uc for the variational parameters and quantities 
such as d and Z increases. 

(3) In the insulating regime, ^'g tends to exhibit gap- 
like behavior in the spin structure factor, 

for small |q|. As seen in the inset of Fig. 14, for t' = 0, 
S{q} is a linear function of |q| at U/t = 0; as U/t in- 
creases, S{q) becomes a quadratic function of |q|, sug- 
gesting that a SC gap opens and becomes large, as will 
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Fig. 13. (Color online) Total energy of projected Fermi sea and 
d-wave states for t' /t = as a function of U/t for several system 
sizes. The Mott critical points for both cases are indicated by 
arrows. Red dashed lines are fitted curves proportional to t/U. 
As expected from Fig. 4(f), the results for \t' /t\ > almost 
coincide with the present curves for t'/t = in the insulating 
regime. 




(0,0) (0,u) q (71,71) (0,0) 

Fig. 14. (Color online) U/t dependence of spin structure factor 
in d-wave state for t'/t = and L = 14 {Uc/t = 6.61). The inset 
shows the magnification of the region near the F point on the 
F-X line, and the symbols are common to the main panel. S(q) 
on the P-M line (nodal direction) exhibits basically the same 
behavior for small |q|. Data points of different Us for U > Uc 
almost overlap one another. 
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Table I. Optimized variational parameters of two wave functions, 

and ^J^, studied in §4, in respective insulating regimes (U > 
Uc) for various t' /t. The final digits may include some errors. 
L = 14. 





\t'/t\ 


9 


Ad/t 




m' 


C/t 




0.0 


0.266 


1.291 


0.831 




0.0 


(U = 


0.25 


0.254 


1.339 


0.840 


-0.0677 


0.200 


7.5t) 


0.3 


0.254 


1.335 


0.837 


-0.0665 


0.230 




0.4 


0.254 


1.349 


0.837 


-0.0690 


0.279 




0.5 


0.253 


1.355 


0.832 


-0.0688 


0.317 




0.0 


0.140 




0.923 






{U = 


0.25 


0.119 




0.904 


0.0776 




12t) 


0.4 


0.121 




0.829 


0.1317 





Table II. Three energy components, total energy and spin struc- 
ture factor at AF wave number, calculated using and in 
respective insulating regimes {U > Uc) for various t' /t. L = 14. 
The final digits may include some errors. Because each energy 
component and E/t are averaged independently, E/t does not 
precisely coincide with the sum of its components. 



I' 


\t'/t\ 


Et/t 


Ef/t 


Eij/t 


E/t 


S(G) 




0.0 


-0.6922 


0.0 


0.2548 


-0.4384 


7.54 


{U = 


0.25 


-0.6888 


-0.0008 


0.2533 


-0.4375 


7.28 


7.5t} 


0.3 


-0.6899 


-0.0012 


0.2532 


-0.4370 


7.27 




0.4 


-0.6848 


-0.0019 


0.2529 


-0.4357 


7.06 




0.5 


-0.6829 


-0.0033 


0.2522 


-0.4340 


6.90 


^'<^ 
{U = 
12t) 


0.0 
0.25 
0.4 


-0.4383 
-0.3455 
-0.3276 


0.0 
-0.0018 
-0.0053 


0.1809 
0.1576 
0.1573 


-0.2575 
-0.1899 
-0.1752 


15.93 
3.81 
2.09 



be mentioned in §5. For U > Uc, the quadratic behavior 
of S'(q) becomes clearer; it is possible that a spin gap 
opens in the insulating regime. As discussed later, the 
frustration makes no difference at this point. This gap- 
like behavior is in contrast to the case of , as will 
be argued in §4. Strictly, however, the insulating state 
represented by can be gapless in the spin sector, in 
the same manner that the d-wave SC is gapless in the 
node direction. To settle this point, further studies are 
necessary. 

(4) In the preceding study for the anisotropic trian- 
gular lattice,^ a band rcnormalization effecf*^ is taken 
into account by optimizing t' in as a variational pa- 
rameter, independently of t' given in the Hamiltonian. In 
the insulating regime, the effective t' is significantly re- 
duced to an almost nonfrustrated value, namely t' /t ~ 0. 
Thereby, the Fermi surface almost recovers the nesting 
condition for the square lattice, leading to highly devel- 
oped short-range AF correlation. In the trial states stud- 
ied here, the band rcnormalization effect is not included, 
but the various results are quantitatively similar to those 
of the previous study, and the AF correlation consider- 
ably increases for U > Uc, as shown in Figs. 14 and 25. 
This result is mainly caused by the behavior of the chem- 
ical potential Q [Fig. 4(e)], which changes its sign for 
U > Uc so as to recover the nesting condition, instead of 
by the band rcnormalization. 

(5) In Table I, we show the optimized parameters in 
the insulating regime {U > Uc) for various t'/t. Note 



that the parameters vary only very slightly, except for 
(/t; the optimized \1/q is almost unchanged with vary- 
ing t'/t. In Table II, we list the total energy and energy 
components for U > Uc, calculated using the optimized 
^Q. E is again almost independent of t'/t, because Ef 
makes a very slight contribution, even for large t'/t (see 
also Table III). As shown in the final column of Table 
II, the AF spin correlation retains a considerably large 
value for large t'/t. This indicates that is stabilized 
by preserving the nesting condition for the square lattice; 
in other words, the quasi-Fermi surface is retained at the 
gap maxima (tt, 0) and (0, tt), at the cost of the energy re- 
duction due to the diagonal hopping or frustration. Thus, 
the AF correlation is a key factor for stabilizing ^q. 

Regarding the AF correlation, it is known for the SC 
states with rf-type symmetries that the gap maxima over- 
lap with the hot spot, namely, the intersection of the 
Fermi surface and the magnetic Brillouin zone bound- 
ary. Thus, it is possible that the shape of the gap func- 
tion Ak deviates from that of the simple d-wave, particu- 
larly for large t'/t. This is an interesting future problem. 

4. Mott Transitions in Projected Fermi Sea 

In this section, we discuss the Mott transition in 

as a continuation of the previous studies for t'/t = 0.^'"^* 
This transition has features different from those of "^q, 
although always has a higher energy than \1/q within 
the present model. In §4.1, we make a careful analysis 
for t' = 0. In §4.2, we consider the t'/t dependence, and 
contrast the properties in with those in '^q. 

4.1 Case for t' = 

The existence of a Mott transition in "i/^ for t' /t = 
was first pointed out in ref. 38, in which the critical value 
was estimated as Uc/t ^ 8.8 by analyzing the cusplike 
behavior in energy components and the discontinuity in 
n(k) for L = 10 and 12. For these sizes, the transition 
appeared continuous. Certainly, we can find cusplike be- 
havior in Fig. 13 [or in the main panel of Fig. 15(d)], 
where the U /t dependence of total energy is shown. Be- 
cause the system-size dependence is considerably large 
near the cusp, we show the magnification in the vicin- 
ity of the cusp in the inset of Fig. 15(d). For L < 14, 
E/t is a smooth function of U/t and has a unique op- 
timized value, whereas for L = 16 and 18, E/t exhibits 
hysteresis or double-minimum behavior near the critical 
point, Uc/t = 8.59 and 8.73, respectively, as in the case 
of \E'q studied in the preceding section. Correspondingly, 
the two variational parameters, g and /U, exhibit discon- 
tinuities at Uc for L > 16, as shown in Figs. 15(a) and 
15(b), although their magnitudes are small compared 
with those of *q in Figs. 4(a) and 4(c). Thus, this tran- 
sition, at least for t'/t = 0, is ascertained to be of the 
first order. 

In Fig. 16, the doublon density is plotted as a func- 
tion of U/t. At the critical point Uc, the order parameter 
d of Mott transitions has a discontinuity for L > 16, 
and suddenly drops to a small value. Simultaneously, the 
doublon-holon binding parameter /i becomes large and 
approaches 1, as shown in Fig. 15(b). These results again 




Fig. 15. (Color online) (a)-(c) Optimized variational parameters (at global minima of E/t) for projected Fermi sea as a function of 
U/t: (a) Onsite correlation (Gutzwiller) factor, (b) Doublon-holon binding factor between nearest-neighbor sites. The insets in (a) and 

(b) are magnifications near Uc for t' /t = 0. The data of local (but not global) minima for L = 16 and 18 are added using open symbols. 

(c) Same factor between diagonal-neighbor sites, (d) Optimized total energies for t' /t = 0, —0.25 and —0.4. Data for several system 
sizes are compared. The inset represents a magnification near Uc for t' /t = 0. For the data of L = 16 and 18, open (closed) symbols 
are used for the metallic (insulating) regime to emphasize the hysteresis. The critical values indicated by arrows are Uc/t = 8.59 and 
8.73 for L = 16 and 18, respectively. The other symbols and the abscissa scales are common to all the main panels. In (c), we add the 
data of including fi' for t' /t = 0. Although fi' has a finite value in the metallic region, its effect is almost negligible. 



suggest that this transition is a Mott transition with the 
doublon-holon binding mechanism similar to that of '^q. 




U/t 

Fig. 16. (Color online) Behavior of doublon density d near the 
Mott critical points in the projected Fermi sea "Jq^ for three 
values of t' /t. Data of several system sizes are plotted for each 
t' /t. For t' /t = and L = 16 and 18, the critical values of 
first-order transitions are indicated by arrows. 



Actually, we have verified that the feature of electron 
configurations is considerably different for U < Uc and 
for U > Uc (not shown), as in Fig. 7. 

The transition is corroborated by the behavior of n(k) 
and iV(q), in the same way as for the c?-wave state. In 
Fig. 17, n(k) [eq. (16)] is plotted along the path F-X-M- 
F for various values of U/t. Because is metallic for 
U < Uc, the Fermi surface can be defined in any direc- 
tion; in this path, kp points are located at the X point 
and at the midpoint of the F-M segment. For U < Uc, 
the discontinuities are apparent at both points of kp, 
whereas for U > Uc, the discontinuities vanish, indicat- 
ing a gap opening. In Fig. 18, we show the magnitude 
of discontinuity, Z, measured at the X point, because 
the extrapolation error is small. As U/t increases, the 
quasi-particle renormalization factor, Z, monotonically 
decreases, and vanishes at [4; the effective mass diverges, 
and ^fg^ becomes insulating. For L < 14, Z is continuous 
near Uc, as reported in Fig. 17 in ref. 38, whereas Z has 
an appreciable discontinuity at Uc for L > 16 and prob- 
ably for L = oo. Incidentally, the discontinuity in Z in 
is smaller than that in (Fig. 10); the character 
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of the first-order transition is less conspicuous in ^'g^. In 
contrast, in dynamical mean- field theory for the hyper- 
cubic lattice, Z continuously decreases and vanishes at 
U — Uc without a discontinuity. 

In Fig. 19, A^(q) [eq. (17)] is plotted for various values 
of U/t. The behavior for small |q| is basically the same 
as that for *q [Fig. 11], namely, A^(q) oc |q| for U < Uc, 
whereas N{q) tends to be proportional to q^ for U > Uc- 
This follows that a charge gap opens for U > Uc- 

From all the above results, we regard this transition 
as a Mott (metal-to-nonmagnetic-insulator) transition, 
which is caused by the doublon-holon binding mecha- 
nism, basically in the same way as in '^q- 

Finally, we consider the spin degree of freedom. In 
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Fig. 17. (Color online) U/t dependence of momentum distribu- 
tion function of for t' /t = 0. In this case, Z (Fig. 18) is 
measured at kp indicated by an arrow on the upper axis (X 
point). 




U/t 

Fig. 18. (Color online) Quasi-particle renormalization factor Z 
of projected Fermi sea, estimated from discontinuities of n(k) 
on X-M line in Brillouin zone. Data for three values of t' /t and 
various system sizes are plotted as a function of U/t. The tails 
for U > Uc are mainly caused by finite-sized effects. 



Fig. 20, S'(q) [eq. (18)] is plotted for various values of 
U/t. As U/t increases, the AF correlation S{G) increases, 
and increases abruptly near the critical point. Although 
the behavior of S{q) is, as a whole, similar to that of 
(Fig. 14), S{G) is twice as large for ^'g^. However, the 
sublattice magnetization. 



1 



(19) 



remains zero within the statistical fluctuation. Thus, a 
long-range order is not realized, although the AF corre- 
lation is considerably enhanced in the insulating regime. 
Note that the behavior of 5(q) for |q| is linear in 
|q| for arbitrary U/t, as shown in the insets of Fig. 20, 
in contrast to the case of 5'g (inset of Fig. 14). This 
strongly suggests that the spin gap is absent. Thus, 
for U > Uc represents a nonmagnetic insulator with- 
out a spin gap, which is considered to be realized in k- 
(ET)2Cu2(CN): 



58 



4-2 Effect of frustration 

First, we consider the effect of frustration on the prop- 
erties of the transition. In Figs. 15(a)-15(c), we show 




Fig. 19. (Color online) U/t dependence of charge density struc- 
ture factor of for t' /t = 0. The symbols denote the same 
values of U/t as in Fig. 20. Uc/t = 8.59. 



3 

Co 




Fig. 20. (Color online) U/t dependence of spin structure factor 
of projected Fermi sea for t' /t = 0. Uc/t = 8.59. The insets are 
the close-ups near the F point on (a) the F-M line and (b) the 
F-X line. 
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Table III. Comparison of ratio p = E^i/Et among the different 
phases indicated in the brackets for the d-wave state and the 
projected Fermi sea. The abbreviations 'ins.' and 'wSC denote 
insulating and weak SC, respectively. The systems of L = 14 is 
used. Values of p are shown as percentages. 





\t'/t\ 


U/t 







6.25 


7.5 


13 




0.25 


3.3 


1.4 (SC) 


0.12 (ins.) 






0.4 


8.8 


7.6 (wSC) 


0.28 (ins.) 






0.25 


3.3 


3.3 (metal) 


3.2 (metal) 


0.4 (ins.) 


0.4 


8.8 


8.8 (metal) 


8.5 (metal) 


1.1 (ins.) 



the optimized variational parameters for t' /t = —0.25 
and —0.4, as well as t' — 0, for systems up to L = 16. 
Let us consider the doublon-holon binding parameter /i 
[Fig. 15(b)] as a typical parameter. For small systems 
(e.g., L — 10), /i is a smoother 'S'-shaped function of 
U/t than that for t' — 0. As the system size increases, 
/i abruptly exhibits semicritical behavior at a slightly 
larger U/t than that for t' ^ [Ujt - 10.95 (11.45) 
for t' /t = -0.25 (-0.4) for L = 16]. In contrast to the 
case of t' ~ 0, the cases of t' /t = —0.25 and —0.4 do 
not indicate a first-order-transition-like discontinuity or 
hysteresis, even for L = 16. In comparing the results 
among the three values of i'/i, we notice that the critical 
behavior tends to be more continuous as the frustration 
becomes strong. We do not consider larger systems in 
this work, because the statistical fluctuation around the 
critical points increases rapidly as L increases. However, 
we assume that this size-dependent critical behavior is a 
sign of a first-order transition. 

These features can be seen in the other variational pa- 
rameters [Figs. 15(a) and 15(c)], total energy [Fig. 15(d)] 
doublon density (Fig. 16) and quasi-particle renormaliza- 
tion factor (Fig. 18). The feature that the critical prop- 
erties tend to become continuous as t' /t increases is op- 
posite to that of studied in §3, but similar to that of 
the path-integral-renormalization-group approach.^'' 




Fig. 21. (Color online) U/t dependence of the spin structure fac- 
tor of the projected Fermi sea for (a) t' /t = —0.25 and (b) 
t'/£ = —0.4. The scale of the vertical axis is common to both 
panels. 



Next, we consider the effect of frustration on the insu- 



lating state. As shown in Table I, when t' /t varies, the 
optimized parameters, namely, the wave function, appre- 
ciably changes, in contrast to that of ^'q. Accordingly, 
the physical quantities with respect to ^'q^ vary with 
i'/t, as shown in Table II. For energy, a notable point is 
that the contribution of Et' is still strongly suppressed 
in similarly to that for ^q. To consider quantita- 

tively how Eti behaves when U /t varies, we list the ratio 
p = Et' /Et in the different phases in Table III. p almost 
maintains its value of the noninteracting case {U ~ Q) in 
the metallic phase, but in the insulating phase, p drops to 
a very small value, although not as small as in This 
feature is affected by p'. As seen in Fig. 15(c), p' abruptly 
decreases at Uc, indicating that the doublon-holon bind- 
ing in the diagonal direction, which assists local diagonal 
hopping for large U/t, also becomes less advantageous for 
U > Uc in ^Q^- The decrease in p' is anticipated to en- 
hance AF correlation, mentioned below. 

Returning to Table II, we find that the AF spin corre- 
lation ^(G) markedly increases for t' /t — 0, but abruptly 
decreases as \t' /t\ increases, in sharp contrast with that 
of To elucidate the situation, we plot S{q) for 

t'/t = -0.25 and -0.4 in Figs. 21(a) and 21(b), respec- 
tively (cf. also Fig. 20 for t' = 0). Although the magni- 
tude of S{G) abruptly decreases as \t' /t\ increases from 
to 0.25, the peak position of S'(q) remains at q = G. 
For t'/t — —0.4, however, in addition to the successive 
decrease in magnitude, the peak of S{q) moves to incom- 
mensurate wave numbers near the M point. 

Thus, the effect of frustration is explicitly reflected in 
physical quantities estimated in the insulating state of 

5. da,2_j,2-wave Superconductivity 

In this section, we study the properties of SC arising in 
the d-wave singlet state ^tg. In §5.1, we deduce the area 
where SC appears in the t'-U plane by distinguishing, 
in the condensation energy, between the contributions 
from a SC gap and from an insulating gap. In §5.2, we 
confirm the appearance of SC by directly observing the 
d-wave pairing correlation ftmction. In §5.3, we consider 
the origin of this SC. 

5.1 Condensation energy 

First, to determine the stability of the d-wave singlet 
state we consider its condensation energy given by 

AE^E{^'}^)-E{^1), (20) 

where £■(*) denotes the optimized variational energy per 
site with respect to VP. In Fig. 22(a), AE/t for various 
t'/t is plotted as a function of U/t; in Fig. 22(b), the 
region near the Mott critical points [Uc/t) is magnified. 

Because the behavior of AE depends on the value 
of t'/t, we first consider the weakly frustrated cases 
{t'/t ^ 0.3). As pointed out in ref. 8, for small values 
of U/t ( 5), AE/t is extremely small; at intermediate 
values of U {— C/onsct ~ 5t-6t), AE/t starts to increase 
slowly at first, then abruptly at the Mott critical value U^ 
of where, in some cases, we can observe a mild cusp 
in Fig. 22(b). As U/t increases further, AE/t has a max- 
imum and then slowly decreases. We are now aware (§3) 
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that, for U > Uc, ^'g becomes an insulating state. Hence, 
the marked increase in AE in this regime is considered 
to originate from the insulating c?-wave gap.^^ Conse- 
quently, the region where substantial energy reduction 
occurs owing to a SC gap is restricted to J/onsot < Uc- 
This idea is supported by the behavior of the d-wave gap 
parameter A/t, which exhibits an appreciable increase 
for the corresponding values of U/t and t' /t, as shown in 
Fig. 4(b). Incidentally, the doublon-holon binding param- 
eter fi increases very similarly to A/t [Fig. 4(c)], suggest- 
ing that this binding plays an active role in the c?-wave 
pairing in the nearest-neighbor sites. 

Next, we proceed to the strongly frustrated cases 
{t' /t ^ 0.4). A major feature of AE, different from that 
in the weakly frustrated case, is that there is no substan- 
tial increase for J/onsct ^ U < Uc- As seen in Fig. 22(b), 
we cannot determine C/onset for t' /t = —0.4 and —0.5, ex- 
cept for a special case, t'/t = —0.4 and L = 10.^^ Corre- 
spondingly, the increase in the d-wave gap A/t [Fig. 4(b)] 
is firmly suppressed in the conductive region (U < Uc), 
compared with those in the weakly frustrated cases. The 
behavior of /i [Fig. 4(c)] again follows that of A/t. Thus, 
SC, if there is any, is expected to be weak in this regime. 

Note that the value of U/t at the maximum of AE/t 
in the insulating regime approximately corresponds to 
the Mott critical point of ^q^, Uc^/t, considered in §4. 
Uc^/t also corresponds to the crossover value, at which 
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Fig. 22. (Color online) Condensation energy, AE/t, of d-wave 
singlet state as function of U/t for various values of t'/t and L. 
(a) AE/t over a wide range of U/t. (b) A close-up of AE/t near 
the Mott critical points; the range is indicated by a dashed-line 
box in (a). The symbols are common to (a) and (b). 



the character of SC changes from the interaction-energy 
origin to the kinetic-energy origin.* 
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Fig. 23. (Color online) Real-space pair correlation function of d- 
wave symmetry along path of r, (0,0)-(0,L/2)-(L/2,L/2)-(0,0), 
for various values of U/t. (a) t'/t = and (b) t'/t = —0.4, 
both for L = 14. The scale of the ordinate axis is common to 
both panels. The data are obtained, using VMC calculations for 
U 0, and from the analytic formula for (7 = 0. 



5. 2 Pair correlation function 

To directly confirm the appearance of d^2_y2-wave 
SC, the d-wave nearest-neighbor pair correlation func- 
tion Pd(r) is convenient for the present approach:^" 

1 



r,r'— x,y 



-5(t,t') 



(At(R,)A,,(R, + r)) 



(21) 



where x and y denote the lattice vectors in the x- and 
y-directions, respectively, and Aj.(Ri) is the creation op- 
erator of a nearest-neighbor singlet. 



(22) 



If Pd(r) has a finite value for |r| oo, off-diagonal long- 
range order exists. For finite systems, however, we have 
to appropriately determine long-distance values of P(i(r), 
particularly, in the cases of small U/t, where the corre- 
lation length is long. In Fig. 23, Pd(r) is plotted for two 
values of t'/t. Although Pd{r) for large |r| should van- 
ish for U/t = 0, spikes of sizable magnitude appear near 
r = {L/2, L/2) [(0,L/2)] for t'/t = [-0.4] for this sys- 
tem size.^^ Furthermore, a trace of this spike structure 
remains up to fairly large values of U/t. Thus, for small 
U/t, we should choose Pd(r) which does not have such 
peculiar finite-sized effects. Fortunately, we found that, 
in the noninteracting cases, the magnitude of |Pd(r)| for 
r = (L/2— 1, L/2), which is almost the farthermost point, 
is very small (less than 10"'*) for arbitrary values of t'/t 
and L. Hence, we employ Pd[{L/2 — 1, L/2)] as the large- 
|r| value, P^, for small U/t, namely, < U < Umax, 
with [/max being the value at which Pd(r) becomes max- 
imum. For the strong-correlation regime ([/ > C/max), 
Pd(r) becomes almost constant for |r| > 3,^^'^^ as shown 
in Fig. 23. Hence, in this regime, we adopt the average 
of Pd(r) for |r| > 3 as P^°°. 
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Fig. 24. (Color online) Large-|r| value of d-wave pair correlation function Pd(r) for (a) small \t' lt\ (< 0.25) and (b) large \t' lt\ (> 0.3). 
Data of various system sizes for each value oit' /t are plotted, (c) Magnification of PJ^ near the Mott critical points. The data for only 
t' /t = 0, —0.25 and —0.4 are shown for clarity. The range of (c) is indicated by a dashed-line boxes in (a) and (b). 



In Figs. 24(a) and 24(b), thus obtained is plotted 
as a function of U /t. In the weakly correlated regime 
(f//t ^ 4), the increase in is small, irrespective of the 
value of t' /t. For weakly frustrated cases {\t' /t\ ^ 0.3), 
P^ starts to increases appreciably at U ^ t^onsot as U/t 
increases, has a peak at U/t — 6-6.25, and then abruptly 
decreases at the Mott critical point U = Uc- The system- 
size dependence of P^ near the peak is weak for L > 12. 
Thus, in these cases, robust d-wave SC certainly occurs 
for ?7onsct ^ U < Uc- On the other hand, for strongly 
frustrated cases [\t'/t\ ^ 0.4 in Fig. 24(b)], no sizable 
increase in P^ is observed at the value corresponding to 
^^onsot- Pd° slowly and monotonically increases until it 
abruptly drops at C/ = Uc-^^ Moreover, the system-size 
dependence is very large. Eventually, robust d-wave SC 
occurs in a limited area, C/onsot < Uc and \t' /t\ ^ 0.3, 
within ^'g. A similar result has been recently obtained 
using a fluctuation exchange approximation.^'* 

In Fig. 24(c), we show the magnification of P^ near 
the Mott critical points. In the insulating regime (C/ > 
Uc), P^ becomes almost independent of t'/t, as men- 
tioned in item (5) in §3.3, decreases rapidly as the sys- 
tem size L increases, and probably vanishes in the limit of 
L — !■ oo. Because the statistical fluctuation in the VMC 
data is much smaller in the insulating regime than in the 
conductive regime, the data are more reliable. The disap- 
pearance of P^ for U > Uc is expected in an insulating 
state. 

5. 3 Properties of superconductivity 

First, we study the relation between d-wave SC and 
AF correlation. In Figs. 14, 25(a) and 25(b), the U/t 
dependence of S'(q) in '^q is shown for t'/t — 0, —0.3 
and —0.4, respectively. As mentioned, robust SC arises 
for t'/t = and -0.3, but does not for t'/t = -0.4 
[Fig. 24]. This is supported by the small-|q| behavior of 
S'(q), shown in the insets of Fig. 25. For t'/t = —0.3, 
S{q) for small |q| tends to be quadratic in |q| as U/t 
increases, indicating that a SC gap develops, whereas for 
t'/t = —0.4, 5(q) remains almost linear in |q|. 

We now focus on the AF wave number G. For U — 0, 



S{q) has a pointed peak for = 0, a rounded peak for 
t'/t = —0.3 and a flat top for —0.4, due to the frustration. 
For the cases oi t' /t — and —0.3, in which robust SC 
appears, 5(G) steadily increases as U/t increases, even 
in the conductive regime, U < Uc [Figs. 14 and 25(a)]. 
In the strongly frustrated case {t'/t — —0.4), in which 
the SC correlation does not develop, the magnitude of 




(0,0) (0,71) q (71,71) (0,0) 

Fig. 25. (Color online) Spin structure factor for (a) t'/t = —0.3 
[Uc/t = 6.69] and (b) t'/t = -0.4 [Uc/t = 6.78] for various U 
near Uc in ^g. The system size is L = 14. The insets show the 
magnification of the region near the F point on the F-X line, 
and the symbols are common to the main panel. Data points of 
different Us for U > Uc almost overlap one another. Similar data 
for t' /t = are given in Fig. 14. 
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S{G) reaches no more than half that for t'/t ~ —0.3, al- 
though S'(q) increases slightly for U < Uc [Fig. 25(b)]. To 
consider the t'/t dependence of S'(q) explicitly, we plot 
S{q) for various values of t'/t in Fig. 26(a) at U/t — 6, 
near which has a maximum (see Fig. 24). When 
t'/t is varied, the change in S'(q) is quantitatively in- 
significant, except near the M point. As \t' /t\ increases, 
S{G) sharply decreases, particularly at t'/t ~ —0.3, 
and G is no longer a characteristic wave number for 
\t'/t\ ;> - 0.45. In Fig. 26(b), we compare the t'/t de- 
pendence of S{G) with that of P^. In respective system 
sizes, when S{G) abruptly decreases for \t' /t\ > 0.25, 
P^ similarly decreases. In Fig. 27, we show the U/t de- 
pendence of <S'(G). Although in every case S{G) gener- 
ally increases as U/t increases, the range of significant 
increase [>S'(G) ^ 2] in the conductive regime roughly 
corresponds to f/onsct ^ U < Uc, and is accompanied by 
a marked increase in P^ (Fig. 24). We have confirmed, 
for a wide range of model parameters, that whenever 
P^ is appreciably enhanced, S{q) has an evident peak 
at q = G. This result strongly supports the idea that 
the SC in this model is induced by AF spin correlation. 
Incidentally, this mechanism is reflected in the ratio of 
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Fig. 26. (Color online) (a) t'/t dependence of spin structure fac- 
tor for U/t = 6, where PJ" is considerably enhanced particularly 
for small t'/t. (b) Spin structure factor at the AF wave num- 
ber G for C//t = 6 as a function of t'/t is denoted by circles 
(left axis). Simultaneously, the d-wave pair correlation function 
is plotted with squares (right axis). 
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energy components. As shown in Table III {U/t — 6.25), 
when SC is weak {\t' /t\ = 0.4), p is only slightly smaller 
than the noninteracting value, whereas for robust SC 
{\t'/t\ — 0.25), p becomes less than half the value for 
U = 0; the diagonal hopping is considerably suppressed. 
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Fig. 27. (Color online) V /t dependence of spin structure factor 
at AF wave number G for various values oit'/t and L. 



Fig. 28. (Color online) (a) Interaction and (b) kinetic compo- 
nents of condensation energy due to iI/q as a function of in- 
teraction strength for four values of t'/t. Various system sizes 
are plotted. The symbols and scales are common to both panels. 



Next, we consider the energy gain in the SC transi- 
tion. The components of condensation energy are shown 
in Fig. 28, that is, the differences in kinetic and interac- 
tion energies between and the actual expression 
is given in the figure. Here, i?kin = (''^kin) = Et + Et', and 
AEi,in + AEij = AE (> 0) [eq. (20)]. In the SC regime 
([/ < Uc), AE'int (Ai?kin) is always positive (negative), 
regardless of the value of t'/t. This indicates that the SC 
transition is induced by the gain in the interaction energy 
at the expense of kinetic energy. This feature smoothly 
continues to the weak-coupling limit (U/t — > 0), and 
is common to conventional BCS superconductors. Al- 
though the component of energy gain switches to the 
kinetic energy at U = 8t-llt, which broadly corresponds 
to U^^, SC is excluded for U > Uc- The kinetic-energy- 
driven SC is not realized at half filling, in contrast to 
that in the doped cases. * 

Incidentally, this behavior of Aitkin and AEjj is not 
restricted to the d-wave state, but is also found in some 
order-disorder transitions. As an example, we plot, in 
Fig. 29(b), AEi^in (= AEt) and AEu for the AF state 
VpQ^ [eq. (14)], which will be discussed in the next sec- 
tion, for t'/t — 0. The interaction part AEjj makes a 
positive contribution to AE for U/t ^ 8, whereas the 
kinetic part AEt contributes for U/t ^ 7. Hence, the 
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behavior is qualitatively identical to that of '^q. Such 
behavior is also observed in SC and CDW states for the 
two-dimensional attractive Hubbard model. 

6. Phase Diagram and Antiferromagnetism 

Up to this point, we have not considered the competi- 
tion with the AF state, but it is a crucial problem when 
drawing a phase diagram. We thus compare the stabil- 
ity between and the AF-ordered state [eq. (14)] 
into which we do not introduce band renormalization pa- 
rameters to equalize the condition with Whenever 
Aaf is finite at half filling, is insulating and has 

an AF long-range order; the sublattice magnetization m 
[eq. (19)] behaves similarly to Aaf- 

In Fig. 30, we show a phase diagram constructed as fol- 
lows. The boundary between AF-I and P-M (\t' /t] < 0.2) 
is determined by the points where the extrapolated val- 
ues of Aaf vanish. The boundary between AF-I and SC 
(or P-I) is determined by the comparison of the total en- 
ergies. The boundary between P-I and P-M is determined 
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Fig. 29. (Color online) (a) Comparison of the condensation en- 
ergy for t'/t = among (indicated by Q/d), (Q/AF) 
and •PG*d (G/d), for which AE = E(Vg'S>fs) - S(7'G*d)- Sev- 
eral system sizes are plotted to consider the L dependence, (b) 
Kinetic (closed symbols) and interaction (open symbols) parts 
of the condensation energy due to VPg^ for t'/t = 0. Four system 
sizes are used. 
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Fig. 30. (Color online) Ground-state phase diagram in t' -V space 
constructed from VMC calculations performed in this study. 
The abbreviations AF-I, P-I (P-M) and SC denote AF insula- 
tor, paramagnetic insulator (metal) and superconductor, respec- 
tively. 



by C/c/i for the largest systems for each t' jt obtained in 
§3. The last boundary is extended to t' /t = Q if is 
not allowed, as indicated by a dashed line. It is unneces- 
sary to fix the boundary between SC and P-M, because 
the small magnitude of A often survives for extremely 
small U/t] instead, we show J/onsot/^ with a dotted line. 

For t'/t = 0, a continuous metal-to- AF-insulator tran- 
sition occurs at U = U^^ = O,^'' owing to the complete 
nesting condition. This AF state is very stable and con- 
tinues to the Heisenberg limit {U/t ^ oo). Our result is 
consistent with it, as shown in Fig. 29(a), where the con- 
densation energy for is always larger than that for 
The nature of the continuous transition seems to be 
preserved in the region of small \t'/tl although f/^F be- 
comes finite, and the tendency toward a first-order tran- 
sition gradually develops as \t' /t\ increases. Note that as 
the frustration becomes strong, \1/qF is rapidly destabi- 
lized, and surrenders to ^'g at t' = t'^ ^ —Q.2t. More- 
over, \t'^/t\ tends to decrease as U/t increases; this fea- 
ture is consistent with the result that the AF state has 
the largest range of n at t7 ^ W \n a, phase diagram in 
the n-U plane. * However, this tendency is not in accord 
with the arguments of the J -J' spin model, which is an 
effective model of eq. (1) for large U/t. Various studies of 
the J -J' model^^ concluded that the AF order vanishes 
at much larger values: J'/J - 0.4 {\t' /t\ ~ 0.63). 

We consider that this discrepancy is primarily due to 
the choice of variational states: (1) ^'g does not have 
a seed of AF long-range order, although the AF short- 
range correlation appreciably develops, and (2) in ^'g^, 
we have not allowed for the band renormalization effect. 
These two requirements are satisfied simultaneously by 
adopting a coexisting state of AF and SC gaps,^'''^^ with 
a band renormalization effect. In fact, we have per- 
formed VMC calculations using such a wave function for 
the anisotropic triangular lattice and found that the area 
of the AF phase considerably expands. A similar result 
has been obtained independently by Chen,^° and also for 
the checkerboard lattice by Koga et al.^^ Thus, it is ur- 
gent that the wave function is improved in this line to 
refine the phase diagram. 

7. Conclusions 

7. 1 Summary 

Using an optimization variational Monte Carlo 
method, we have studied the half-fiUed-band Hubbard 
model on frustrated square lattices, given by eq. (1). 
Our primary aim is to understand the mechanisms of 
the Mott transition and of the _y2-wave SC arising in 
the Hubbard model. To this end, we introduce an inter- 
site correlation factor that controls the binding between 
a doublon and a holon into the trial functions: normal 
(Fermi sea), d-wave singlet and AF states. We have suc- 
ceeded in describing the d-wave SC and a Mott transition 
simultaneously in a single approach. We itemize our main 
findings: 

(1) Within the d-wave singlet state, a first-order Mott 
(conductor-to-nonmagnetic-insulator) transition occurs 
at Ucj which is approximately the bandwidth, for arbi- 
trary t'/t. In the insulating regime, most doublon-holon 
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pairs are actually confined within the nearest-neighbor 
sites, in contrast to the case in the conductive regime. 
The critical Uc/t gradually increases as the frustration 
becomes strong. This transition is not directly related to 
a magnetic order. 

(2) We have confirmed that in the projected Fermi 

first-order Mott transition without relevance to 
magnetism also arises at a larger U/t than the bandwidth 
for arbitrary t' /t, although the state does not have the 
lowest energy. 

(3) The nonmagnetic insulating state (d-wavc singlet 
state for U > Uc) has a considerably low energy and 
a strong short-range AF correlation. According to the 
small-|q| behavior of S'(q), the d-wave state tends to have 
short singlet bonds owing to the nearest-neighbor pairing 
[eq. (12)], in contrast to the projected Fermi sea, which 
clearly does not have a spin gap. 

(4) Robust SC with d^2_y2-wa.ve symmetry ap- 
pears for moderate values of U/t (~ 6) and t'/t 
(0.2 ^ \t'/t\ ^ 0.35). This area is adjacent to both do- 
mains of a Mott insulator and an AF long-range order. 
The phase diagram obtained in this study is shown in 
Fig. 30. 

(5) By comparing the pair correlation function with 
S'(q), it is found that robust SC is always accompanied by 
appreciably enhanced short-range AF spin correlation, 
which is weakened by the frustration and almost vanishes 
for \t'/t\ ^ 0.35. The SC transition is induced by the 
gain in interaction energy; this mechanism is identical 
to that in the weak-correlation limit as well as that of 
conventional BCS superconductors. 

(6) The AF long-range order prevailing in the weakly 
frustrated cases {\t' /t\ ^ 0.2) is rapidly destabilized as 
\t'/t\ increases if a band renormalization effect is not in- 
troduced. 

7.2 Further discussions 

(1) In comparing the present study for the frustrated 
square lattice [Fig. 1(a)] with the preceding study,^ in 
which almost the same wave hmctions are applied to 
the anisotropic triangular lattice [Fig. 1(b)], the results 
for the two lattices are qualitatively identical, indicating 
that the two types of frustration work similarly unless 
\t'/t\ is too large. However, the critical values with re- 
spect to t'/t are approximately doubled for the latter 
lattice; namely, the AF state becomes unstable at ap- 
proximately \t' /t\ = 0.2 for the former and 0.4 for the 
latter, and the robust SC disappears at approximately 
\t'/t\ = 0.35 for the former, and 0.8 for the latter. This 
can be explained by the difference in the number of frus- 
trated bonds. 

(2) Recently, using a VMC method with a two-body 
long-range Jastrow factor for the square lattice {t' = 0), 
Capello et al7^ found that a metal-to-insulator transi- 
tion arises, similar to that in ^'g^, for example, in the 
behavior of Z and d. The critical value of their function 
is Uc/t ~ 8.5, which is close to 8.73 {L = 18) in ^-^s. 
Although their transition is regarded as continuous, evi- 
dence of the first order is possibly be found by a detailed 
analysis of larger systems. In their wave function, no ex- 
plicit (four-body) doublon-holon binding factor is intro- 



duced, but a short-range part of the two-body Jastrow 
factor may substantially work as a binding factor under 
the condition of strong electron repulsion at half filling. 
It will be interesting to reveal the relation between the 
two wave functions. 

(3) In this paper, we have restricted the electron den- 
sity to half filling (n = 1). When carriers arc doped, 
unless the doping rate |1 — n| is too large, the doublon- 
holon-binding effect remains significant, as we showed 
for t' = in the previous letter.^ The Mott transition at 
half filling changes to a crossover from weakly to strongly 
correlated regimes. As [1 — n[ increases, the AF order is 
rapidly destabilized, and the SC phase expands to the re- 
gion of large U/t, consistent with the behavior of high-Tc 
cuprates. We will report a detailed description of doped 
cases with the effect of frustration elsewhere. 
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